Jammed particulate systems are inherently nonharmonic 
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Jammed particulate systems, such as granular media, colloids, and foams, interact via one-sided 
forces that are nonzero only when particles overlap. We find that systems with one-sided repulsive 
interactions possess no linear response regime in the large system limit (N — >• oo) for all pressures 
p (or compressions A0), and for all N near jamming onset p — »> 0. We perform simulations on 2D 
frictionless bidisperse mechanically stable disk packings over a range of packing fractions = (/)—(/) j 
above jamming onset <fij. We apply perturbations with amplitude S to the packings along each eigen- 
direction from the dynamical matrix and determine whether the response of the system evolving 
at constant energy remains in the original eigenmode of the perturbation. For 5 > 8 C , which we 
calculate analytically, a single contact breaks and fluctuations abruptly spread to all harmonic 
modes. As S increases further all discrete harmonic modes disappear into a continuous frequency 
band. We find that (5 C ) ~ A(f)/N x , where 1 > A > 0.5, and thus jammed particulate systems are 



inherently nonharmonic with no linear vibrational response regime as N 
A(/>, and as A0 — >• at any N. 

PACS numbers: 83.80.Fg,63.50.-x,62.30.+d,61.43.-j 
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Introduction Granular materials, which are collec- 
tions of macroscopic grains that interact via contact 
forces, such as sand, powders, pharmaceutical, and 
consumer products, display strongly nonlinear spatio- 
temporal dynamics even when they are weakly driven. In 
stark contrast to conventional solids [l| granular solids 
possess nonamne, hysteretic, and time-dependent me- 
chanical response [2[, and dispersive, attenuated, and 
noisy acoustic response [3|, |4[ for micro-strains. 

Crystalline and amorphous atomic and molecular 
solids display well-defined linear response regimes for 
small perturbations. Similarly, there has been a large 
research effort to identify linear response regimes for 
granular and other jammed particulate systems. Ex- 
amples include effective medium theory [5| for granular 
media, which provides predictions for the elastic mod- 
uli as a function applied pressure, and approaches that 
assume the vibrational modes of static, mechanically sta- 
ble (MS) packings obtained from the dynamical matrix 
in the harmonic approximation describe the mechanical 
response [6j , vibrations [3, [8j , and heat flow [9j of weakly 
perturbed and fluctuating particulate systems. 

However, it has not been determined whether jammed 
particulate systems possess a linear response regime, and 
if so, over what range of perturbation amplitudes and 
timescales. To address this fundamental question for 
granular media, it is important to understand separately 
the manifold contributions to nonharmonicity includ- 
ing nonlinear, dissipative, and frictional particle interac- 
tions [lo|| , inhomogeneous force propagation [lll[l2j , and 
breaking and forming of intergrain contacts |3j- Here, 
we describe computational studies to quantify perhaps 
the most important contribution to nonharmonicity in 



jammed particulate media — the one-sided nature of con- 
tact interactions — interparticle forces are only nonzero 
when two grains are in contact, but are strictly zero when 
they are out of contact. 

We find that one-sided interactions make jammed par- 
ticulate materials inherently nonharmonic, i.e. nonhar- 
monic even in the limit of vanishing perturbation am- 
plitude, due to changes in the contact network following 
the perturbation [13]. Specifically, we employ the har- 
monic approximation and calculate the eigenmodes of the 
dynamical matrix [l4[ for MS frictionless packings, sub- 
ject the packings to vibrations along the harmonic set of 
eigenmodes, and quantify the frequency content of the re- 
sponse versus the perturbation amplitude S. We find that 
systems become nonharmonic (i.e. the response is not 
confined to the original mode of excitation) when only a 
single contact is broken (or gained) at a critical S c that 
depends on the original mode of excitation. For S > S c 
the response first spreads to all (harmonic) eigenmodes 
with an amplitude that scales inversely with frequency, 
and then becomes continuous with an average frequency 
that decreases with 5. We show that (S c ) averaged over 
the modes of excitation tends to zero in the large system 
limit even for highly compressed systems, and tends to 
zero in the limit of zero compression at all system sizes. 
Thus, jammed particulate systems possess no harmonic 
regime in the large system limit and at jamming onset 
for any system size. 

Model and Simulations We focus on frictionless MS 
packings of bidisperse disks in 2D with system sizes in the 
range N = 12 to 1920 particles using periodic boundaries 
in square simulation cells (2N/3 disks with diameter <j 
and 7V/3 disks diameter 1.4<j). The disks interact via the 
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FIG. 1: (a) Mechanically stable (MS) packing of frictionless 
disks for N = 12 at = 10 -5 (black solid) and a packing 
perturbed along the 6th eigenmode of the dynamical matrix 
by S = 0.1a (red dashed). The vector lengths are proportional 
to the displacements, (b) An intensity plot of the logarithm 
of the power spectrum |.R(cj)| 2 as a function of frequency uj 
and perturbation 5 along the 6th eigenmode of the system in 
(a) after 170 oscillations. The solid horizontal lines indicate 
the 22 harmonic eigenfrequencies for (a). The inset shows 
the same calculation except for a two-sided linear spring po- 
tential, (c) Same as (b) except for N = 58 at A0 = 10 -5 
with perturbation in mode 40 after 150 oscillations. The in- 
set shows a close-up of the transition. 

linear repulsive spring potential 

where is the center-to-center separation between disks 
i and j, e is the characteristic energy scale, Q(x) is the 



Heaviside function, and = (<j.j + Oj)/2 is the average 
diameter. We have also studied systems with Hertzian 
and purely repulsive Lennard- Jones interactions, but the 
repulsive linear spring potential provides a 'lower bound' 
on the degree of nonlinearity arising from one-sided in- 
teractions. Energy, length, and timescales are measured 
in units of e, <r, and y^m/ecr, respectively. 

The MS packings were generated using the com- 
pression and energy minimization protocol described in 
Ref. [16]. Each MS packing is characterized by a pack- 
ing fraction above which the potential energy V and 
pressure p of the system begins to increase from zero. 
The distance in packing fraction from (j)j is tuned from 
A(j) = 10 -8 to 10 _1 and the positions of the particles are 
accurate to 10 -16 at each A<p. We calculate the eigen- 
frequencies Ui and eigenmodes ei = {ef , ef , . . . , ef} = 

{ e l^ e l^ e li^ e l^--^ e x^ e yi} ( with = 1) in the har- 
monic approximation from the dynamical matrix evalu- 
ated at the MS packing. Since the systems are mechani- 
cally stable, the j\f — 27V' — 2 eigenfrequencies uji > [Uj], 
where N f = N — N r and N r is the number of rattler parti- 
cles with less than three contacts per particle. We index 
the eigenfrequencies from smallest to largest, i = 1 to jV, 
removing the two trivial eigenfrequencies corresponding 
to uniform translations. 

To test whether the packings possess a harmonic 
regime, we apply displacements to individual particles 
and then evolve the system at constant total energy E. 
Specifically, at time t = 0, we apply the displacement 

R-R° = 5e u (2) 

where the new configuration R = {R\, R2 • • • , Rn} = 
{xi , y\ , . . . , xn , yjy} , and R° is the original MS packing. 
We remove rattlers from the MS packings prior to apply- 
ing the perturbations. A sample perturbation for N = 12 
along the 6th mode is shown in Fig. [T] (a). For t > 0, we 
solve Newton's equations of motion at constant E, and 
measure the particle displacements and number of con- 
tacts as a function of the number of oscillations n for 
perturbations along each mode k. 

Results In Fig. [T] (b), we show the logarithm of the 
power spectrum \R(uj)\ 2 , where R(cj) — f™ T6 dte ZUJt R(t) 
for n = 170 oscillations, where Tq = 2tt/ujq, as an inten- 
sity plot versus the perturbation amplitude S (along the 
6th mode) and uj for the system shown in Fig. [T] (a) with 
linear repulsive spring interactions. This plot demon- 
strates several key features: (1) There is an extremely 
sharp onset of nonharmonicity at log 10 S% /a ~ —6.8. For 
5 < 5%, the system vibrates with uj = ujq. Although 5% 
depends on the excitation mode, the transition for each 
mode is sharp. (2) For 5 > 5%, the response spreads to 
include other harmonic eigenfrequencies (shown as solid 
horizontal lines in Fig. [T] (b)) and \R(uj)\ 2 ~ oj~ 2 simi- 
lar to equipartition in thermal equilibrium. (3) For larger 
perturbations, the power spectrum develops a continuous 
frequency band in which the harmonic eigenfrequencies 
are completely lost. At sufficiently large amplitudes, the 
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FIG. 2: (a) Amplitude A% (solid) and deviation of the number of contacts AiV c = N® — (N c )t relative to the unperturbed 
number N® (dashed) versus perturbation amplitude 5 along four eigenmodes (labeled by mode number) for the system in Fig. [1] 
(c) after n = 10 4 oscillations. The vertical dot-dahsed line indicates 5 = 5% for k — 114. (b) The measured 5%(k) at which 
A k and (N c )t begin to deviate from 1 for perturbations along all eigenmodes at n — 10 4 , N — 58, and A0 = 10 -5 versus the 
calculated deformation amplitude 5 c (k) in Eq. [4] at which the first contact breaks. The inset displays the values 5* — S c (k) at 
which (A%) decays to 0.2 as a function of n. The solid line has slope —1. (c) A 1 ^ versus 5 — S c (k) for each k (open symbols) 
and (A*:) averaged over k (lines) for n — 1 (circles, solid line), 10 2 (squares, dashed line), and 10 4 (triangles, dot-dashed line) 
oscillations after the perturbation. 



dominant contribution to the broad power spectrum ap- 
proaches uj = 0. Note that this crossover to nonharmonic 
frequency response occurs at much larger amplitudes in 
systems with smooth nonlinear interaction potentials. 

For larger systems the transition from harmonic to 
nonharmonic behavior is similar (Fig. [1] (c)). The in- 
set to Fig. [1] (c) shows that large systems display an 
intermediate nonharmonic regime in which a subset of 
harmonic eigenmodes are populated at the onset of non- 
harmonicity, 5 = 5%. To put the effects of one-sided po- 
tentials into perspective, we compare these results with 
those from two-sided spring potentials (i.e. Eq. [I] with 
the argument of replaced by 1 — Rij/crij). For N = 12, 
the transition for systems with one-sided repulsive spring 
interactions occurs at perturbations more than four or- 
ders of magnitude smaller than those for systems with 
double-sided spring potentials [15] and the transition oc- 
curs slowly over a decade in 5 (inset to Fig. [U(b)). 

To quantify the harmonic to nonharmonic transition, 
we calculate the number of particle contacts (N c ) t av- 
eraged over time and define a harmonicity parameter 
A\ that measures the spectral content of the particle 
displacements in the eigenmode direction k at eigenfre- 
quency uok following a perturbation along eigenmode k: 

A k = Jp Tk A ^ft) • g fe cos(uj k t)dt 
5 f™ Th cos 2 (ojkt)dt 

where AR(t) = R(t) - (R(t)) t . A% = 1 for harmonic 
systems and A% « for nonharmonic systems that do 
not oscillate in mode k at uok- We also calculate the har- 
monicity parameter (A%) averaged over all individually 
perturbed modes k. 

In Fig. [2] (a), we plot A\ and the deviation in the time- 
averaged number of contacts AN C = N® — (N c ) t rela- 
tive to the unperturbed value N® versus 5 along several 
modes k for the system in Fig.[T](c). We find that A% for 
each mode k begins to decrease from 1 at the same 5%(k) 
where the average number of contacts (N c ) t begins to 



deviate from N®. For perturbations along each mode fc, 
the transition from harmonic to nonharmonic behavior 
occurs when a single existing contact breaks. To ver- 
ify this, we plot in Fig. [2] (b) 5%(k) versus the predicted 
amplitude 5 c (k) at which the first contact breaks. The 
predicted value 5 c (k) is obtained by solving R? { - = afj for 
all contacting pairs of particles i and j for a given MS 
packing and perturbation along mode fc, and identifying 
the minimum 5 c (k) = min^- where 




(4) 

We find that the 5 at which A% begins to decrease, 
5%(k) = 5 c (k) at which a single contact breaks (as shown 
in Fig. [2] (b)) over a wide range of A</> and N with a rel- 
ative error less than 10 -3 over four orders of magnitude 
in 5 c (k). For larger system sizes, it is possible that new 
contacts can form before existing contacts break, but we 
find that this does not occur for the system sizes and 
compressions studied. 

The rate at which energy input via a perturbation 
along eigenmode k is transferred out of that mode and 
into other displacement modes determines the shape of 
the decay of A%. In Fig. [2] (c), we show A% and (A\) 
versus 5 — 5 c (k) for n = 1, 10 2 , and 10 4 oscillations for 
perturbations along each mode k individually. For small 
n, even though A 1 ^ begins to decrease from 1 at ^ c (^) 7 the 
shape of the decay depends on k and the sharp decrease 
from 1 to occurs at small but finite 5 — 5 c (k). In the 
inset to Fig. [2] (b), we measure the amplitude 5* — 5 c (k) 
at which (A%) decays to a small value (0.2), and find 
5* — 5 c (k) ~ 1/n. For all n and perturbations 5 studied 
in Fig. [2j there is no detectable nonharmonicity from the 
smooth nonlinearities in the potential [15| in Eq. [TJ 

Thus, 5 c (k) is the critical deformation amplitude above 
which MS packings become nonharmonic, and 5 c (k) can 
be calculated exactly using Eq. 2] for each MS packing 
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FIG. 3: (a) Distribution P(S c (k)) (scaled by A<j)/N) versus 8 c {k)N/A<j) for = 10" 2 (dotted), 10" 4 (dot-dashed), and 10" 7 
(dashed) and N = 60 (black), 240 (red), and 1920 (blue). The two solid black lines have slope 1.5 and 2. The inset shows the 
scaling of (6 c (k)) with N for A(j) = 10 -2 (crosses), 10 -3 (pluses), 10 -4 (upward triangles), 10 -5 (diamonds), 10 -6 (squares), 
and 10 -7 (circles). The two solid lines have slope 1 and 0.5. (b) The total energy per particle required to break a single contact 
averaged over k (scaled by A(A(j))(A(j)) 2 ) versus system size N for A<j) = 10 -2 (crosses) 
10 -5 (diamonds), 10 -6 (squares), and 10 -7 (circles). The solid line has slope —1.7. 



10 (pluses), 10 (upward triangles) 



and mode k. In Fig. [3] (a) we show the distribution 
P(S c (k)/a) of critical amplitudes with the vertical and 
horizontal axes scaled by A(j)/N to achieve approximate 
collapse. We find that the distribution P(S c (k)/a) scales 
roughly as a power law (S c (k)/a) a for large 5 C with a « 2 
(1.5) for large (small) A<p. The distribution is cut-off and 
remains nearly constant for S c (k)/cr < A(j)/N. In the in- 
set to Fig. [3] (a), we plot (S c (k)) averaged over k versus N 
over a range of A<p =10 -7 tol0 -2 . As expected from the 
power-law distribution in the main plot, (S c (k)) ~ iV a_1 . 
For all A0, the critical deformation amplitude scales to 
zero in the large system limit. 

The potential energy V of a MS packing prepared at 
A0 is given by V/N = B(A(j)) 2 0, where B is a 0(1) 
constant. We find that the average deformation energy 
E* « ((ukS c (k)) 2 ) for the critical amplitude scales as 
E* ~ A(A^)(A0) 2 /7V /3 , where A(A<j>) is only weakly de- 
pendent on A(j) and ^ « 1.7. For E > E* (labeled A in 



Fig. |3](b)), MS packings are strongly anharmonic. MS 
packings are only harmonic for E < E*, where E* — >• 
in the large system limit for all A<p. Thus, eigenfrequen- 
cies of the dynamical matrix do not describe vibrations 
of MS packings as N — ^ oo for all A<p. 

Conclusion We have shown that one-sided repulsive 
interactions in jammed particulate systems make them 
inherently nonharmonic. In the large system limit at 
any compression and in the A(j) — >• limit at any sys- 
tem size, infinitesimal perturbations will cause them to 
become strongly nonharmonic, which will affect their me- 
chanical response, specific heat, and energy diffusivity. In 
future studies, we will explore the possibility of defining 
dynamic steady-state packings with robust effective har- 
monic modes obtained from average particle positions. 
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